perm filename LPC1.F4[6,ALS] blob
sn#001156 filedate 1972-05-12 generic text, type T, neo UTF8
00100 SUBROUTINE LPC1(IFFY,NPTS,CF,M,RO,ERRN,ERR,SPT,NSP,ISSW)
00200 C THE PARAMETERS IN THE LIST ARE DEFINED AS FOLLOWS:
00300 C IFFY←THE INPUT DATA ARRAY OF INTEGERS...A TIME SERIES
00400 C NPTS←THE NUMBER OF POINTS IN THE TIME SERIES
00500 C CF← THE ARRAY IN WHICH THE FILTER COEFFICIENTS
00600 C ARE TO BE RETURNED.
00700 C M← THE NUMBER OF FILTER COEFFICIENTS DESIRED.
00800 C RO← COMPUTED FILTER SCALE FACTOR.
00900 C ERRN←NORMALIZED ERROR VALUE.
01000 C ERR← ARRAY IN WHICH THE SERIES OF ERROR IS RETURNED.
01100 C SPT← THE ARRAY IN WHICH THE REAL SPECTRUM IS RETURNED.
01200 C NSP← THE NUMBER OF POINTS IN THE REAL SPECTRUM.
01300 C ISSW← AN INTEGER SWITCH USED TO SELECT THE FOLLOWING
01400 C MODES OF OPERATION:
01500 C 1: CALCULATE AND RETURN ALL VALUES IN LIST.
01600 C 2: CALCULATE RO, ERRN, AND CF.
01700 C 3: CALCULATE RO, ERRN, CF, AND ERR.
01800 C 4: CALCULATE RO, ERRN, CF, AND SPT.
01900 C 5: USING CF AND M AS INPUT CALCULATE ERR AND SPT.
02000 C 6: " " " " " " " ERR.
02100 C 7: " " " " " " " SPT.
02200 C 1X: FOR ANY X ABOVE IMPLYING IFFY INPUT,
02300 C SUPPRESS HAMMING WINDOW APPLICATION.
02400 DIMENSION IFFY(1),CF(1),ERR(1),SPT(1),X(512),Y(512)
02500 TPI=2.*355./113.
02600 ISW=ISSW
02700 IHAM=1
02800 C IHAM=1←APPLICATION OF A HAMMING WINDOW.
02900 IF ((ISW -10).GT.0) IHAM=0
03000 IF (IHAM.EQ.0) ISW=ISW-10
03100 C LOC 100 IS CALCULATION OF CF, RO, AND ERRN.
03200 IF (ISW.LE.6) GO TO 10
03300 1 IF (ISW.LT.5) GO TO 100
03400 C LOC 200 IS CALCULATION OF ERR.
03500 2 IF ((ISW.EQ.1).OR.(ISW.EQ.3).OR.(ISW.EQ.5).OR.(ISW.EQ.6)
03600 1 ) GO TO 200
03700 C LOC 300 IS CALCULATION OF SPT.
03800 3 IF ((ISW.EQ.1).OR.(ISW.EQ.4).OR.(ISW.EQ.5).OR.(ISW.EQ.7)
03900 1 ) GO TO 300
04000 RETURN
04100 10 NP=NPTS-1
04200 C CALCULATE THE DIFFERENCE SIGNAL
04300 DO 11 J=1,NP
04400 11 X(J)=IFFY(J+1)-IFFY(J)
04500 IF (IHAM.EQ.0) GO TO 1
04600 C APPLY A HAMMING WINDOW TO THE TIME SERIES.
04700 DT=NP
04800 X(NPTS)=X(NPTS-1)*.08
04900 DO 12 J=1,NP
05000 T=J-1
05100 12 X(J)=X(J)*(.54-.46*COS(TPI*T/DT))
05200 GO TO 1
05300 C CALCULATE AUTOCORRELATION COEFFICIENTS AND SOLVE
05400 C FOR INVERSE FILTER COEFFICIENTS
05410 100 DO 105 J=1,NPTS
05420 105 Y(J)=X(J)
05500 CALL INVFL(Y,NPTS,M,RO,ERRN)
05600 C Y ARRAY RETURNS 1.,A1,A2,...,AM,0.,....0.
05700 MM=M+1
05800 DO 125 J=1,MM
05900 125 CF(J)=Y(J)
06000 GO TO 2
06100 200 MM=M+1
06200 DO 220 J=MM,NPTS
06300 ERR(J)=0.
06400 DO 220 K=0,M
06500 220 ERR(J)=ERR(J)+X(J-K)*CF(K+1)
06600 DO 230 J=1,M
06700 ERR(J)=0.
06800 DO 230 K=0,M
06900 IF ((J-K).GE.1) ERR(J)=ERR(J)+CF(K+1)*X(J-K)
07000 230 CONTINUE
07100 GO TO 3
07200 300 DO 303 J=1,512
07300 X(J)=0.
07400 303 Y(J)=0.
07500 MM=M+1
07600 DO 310 J=1,MM
07700 310 X(J)=CF(J)
07800 C SHUFFLE DATA FOR CALC OF REAL TRANSFORM
07900 MM2=MM/2 +1
08000 DO 320 K=1,NSP
08100 K1=2*K-1
08200 K2=2*K
08300 X(K)=X(K1)
08400 320 Y(K)=X(K2)
08500 NMOD=1
08600 323 IF (2**NMOD.GE.NSP) GO TO 326
08700 NMOD=NMOD+1
08800 GO TO 323
08900 326 IF ((2*MM2).LT.MM) MM2=MM2+1
09000 MOD=1
09100 330 IF (2**MOD.GE.MM2) GO TO 340
09200 MOD=MOD+1
09300 GO TO 330
09400 340 CALL FFTP(X,Y,NMOD,MOD)
09500 CALL RBITS(X,Y,NMOD)
09600 CALL USCRM(X,Y,NMOD)
09700 C CALCULATE LOG OF RECIPROCAL INVERSE FILTER SPECTRUM
09800 DO 350 J=1,NSP
09900 350 X(J)=1./(X(J)*X(J)+Y(J)*Y(J))
10000 DO 360 J=1,NSP
10100 360 SPT(J)=10.*ALOG10(X(J))
10200 RETURN
10300 END
00100 SUBROUTINE INVFL(X,NP,M,RO,ERRN)
00200 DIMENSION X(512)
00300 DIMENSION F(50),TF(50),A(50),R(50)
00400 C CALCULATION OF M+1 LENGTH AUTOCORRELATION SEQUENCE
00500 C FROM THE INVERSE FILTER INPUT SEQUENCE
00600 MP1=M+1
00700 DO 11 JJ=1,MP1
00800 J=JJ-1
00900 MMJ=NP-J
01000 SS=0.
01100 DO 10 I=1,MMJ
01200 IPJ=I+J
01300 10 SS=SS+X(I)*X(IPJ)
01400 11 R(JJ)=SS
01500 RO=R(1)
01600 C FORTRAN IMPLEMENTATION OF LEVINSON'S METHOD FOR
01700 C SOLVING THE AUTOCORRELATION MATRIX
01800 F(1)=1.
01900 ALPHA=R(1)
02000 BETA=R(2)
02100 A(1)=-R(2)/R(1)
02200 GAMMA=A(1)*R(2)
02300 DO 1 N=2,M
02400 NM1=N-1
02500 C=-BETA/ALPHA
02600 IF (N-2) 2,2,3
02700 3 DO 4 J=2,NM1
02800 NN=N-J+1
02900 4 TF(J)=F(J)+C*F(NN)
03000 DO 5 J=2,NM1
03100 5 F(J)=TF(J)
03200 2 F(N)=C*F(1)
03300 ALPHA=ALPHA+C*BETA
03400 BETA=0.
03500 DO 6 J=1,N
03600 NN=N-J+2
03700 6 BETA=BETA+F(J)*R(NN)
03800 Q=-(R(N+1)+GAMMA)/ALPHA
03900 DO 7 J=1,NM1
04000 NN=N-J+1
04100 7 A(J)=A(J)+Q*F(NN)
04200 A(N)=Q*F(1)
04300 GAMMA=0.
04400 DO 8 J=1,N
04500 NN=N-J+2
04600 8 GAMMA=GAMMA+A(J)*R(NN)
04700 1 CONTINUE
04800 C CALCULATE NORMALIZED ERROR ERRN
04900 SM=0.
05000 DO 77 J=1,M
05100 77 SM=SM+A(J)*R(J+1)
05200 ERRN=1.+SM/R(1)
05300 C PLACE COEFFICIENTS IN PROPER LOCATIONS OF THE
05400 C REAL VECTOR X FOR FFT-ING TO GET INVERSE SPECTRUM
05500 DO 51 J=1,M
05600 51 X(J+1)=A(J)
05700 X(1)=1.
05800 MP1=MP1+1
05900 DO 123 J=MP1,512
06000 123 X(J)=0.
06100 RETURN
06200 END
00100 SUBROUTINE FFTP(X,Y,M,L)
00200 DIMENSION X(512),Y(512)
00300 N=2**M
00400 L2=2**L
00500 DO 1 LO=1,M
00600 LMX=2**(M-LO)
00700 LMM=LMX
00800 LIX=2*LMX
00900 SCL=6.283185/LIX
01000 C TEST FOR PRUNING
01100 IF (LO-M+L) 20,30,30
01200 20 LMM=L2
01300 30 DO 1 LM=1,LMM
01400 ARG=(LM-1)*SCL
01500 C=COS(ARG)
01600 S=SIN(ARG)
01700 DO 1 LI=LIX,N,LIX
01800 J1=LI-LIX+LM
01900 J2=J1+LMX
02000 T1=X(J1)-X(J2)
02100 T2=Y(J1)-Y(J2)
02200 X(J1)=X(J1)+X(J2)
02300 Y(J1)=Y(J1)+Y(J2)
02400 X(J2)=C*T1+S*T2
02500 1 Y(J2)=C*T2-S*T1
02600 RETURN
02700 END
00100 SUBROUTINE RBITS(X,Y,M)
00200 DIMENSION X(512),Y(512),L(9)
00300 EQUIVALENCE(L9,L(1)),(L8,L(2)),(L7,L(3)),(L6,L(4)),
00400 1 (L5,L(5)),(L4,L(6)),(L3,L(7)),(L2,L(8)),(L1,L(9))
00500 C PERFORM BIT REVERSAL
00600 DO 70 J=1,9
00700 L(J)=1
00800 IF (J-M) 71,71,70
00900 71 L(J)=2**(M+1-J)
01000 70 CONTINUE
01100 JN=1
01200 DO 60 J1=1,L1
01300 DO 60 J2=J1,L2,L1
01400 DO 60 J3=J2,L3,L2
01500 DO 60 J4=J3,L4,L3
01600 DO 60 J5=J4,L5,L4
01700 DO 60 J6=J5,L6,L5
01800 DO 60 J7=J6,L7,L6
01900 DO 60 J8=J7,L8,L7
02000 DO 60 JR=J8,L9,L8
02100 IF (JN-JR) 61,61,62
02200 61 R=X(JN)
02300 X(JN)=X(JR)
02400 X(JR)=R
02500 FI=Y(JN)
02600 Y(JN)=Y(JR)
02700 Y(JR)=FI
02800 62 JN=JN+1
02900 60 CONTINUE
03000 RETURN
03100 END
00100 SUBROUTINE USCRM(X,Y,M)
00200 DIMENSION X(512),Y(512)
00300 I=2**M
00400 LO2=I/2
00500 SA=3.141593/I
00600 X(1)=2.*(X(1)+Y(1))
00700 Y(1)=0.
00800 LLL=LO2+1
00900 X(LLL)=2.*X(LLL)
01000 Y(LLL)=-2.*Y(LLL)
01100 DO 33 K=2,LO2
01200 J=I-K+2
01300 ANG=SA*(K-1)
01400 C=COS(ANG)
01500 S=SIN(ANG)
01600 AA=X(K)+X(J)
01700 BB=Y(K)-Y(J)
01800 CC=X(K)-X(J)
01900 DD=Y(K)+Y(J)
02000 XR=C*DD-S*CC
02100 XI=S*DD+C*CC
02200 X(K)=AA+XR
02300 Y(K)=BB-XI
02400 X(J)=AA-XR
02500 Y(J)=-BB-XI
02600 33 CONTINUE
02700 RETURN
02800 END